logit = function(prob) {
    log(prob) - log(1 - prob)
}

expit = function(logodds) {
    1/(1 + exp(-logodds))
}


getlogop = function(p0, p1) {
    log(p0) + log(p1) - log(1 - p0) - log(1 - p1)
}

getlogrr = function(p0, p1) {
    log(p1) - log(p0)
}

getatanhrd = function(p0, p1) {
    atanh(p1 - p0)
}


## Function for checking if two things are equal within numerical precision
same = function(x, y, tolerance = .Machine$double.eps^0.5) {
    abs(x - y) < tolerance
}


## Functions for wrapping estimation results into a nice format
WrapResults = function(point.est, cov, param, name, va, vb, converged) {
    
    se.est = sqrt(diag(cov))
    
    conf.lower = point.est + stats::qnorm(0.025) * se.est
    conf.upper = point.est + stats::qnorm(0.975) * se.est
    p.temp = stats::pnorm(point.est/se.est, 0, 1)
    p.value = 2 * pmin(p.temp, 1 - p.temp)
    
    names(point.est) = names(se.est) = rownames(cov) = colnames(cov) = names(conf.lower) = names(conf.upper) = names(p.value) = name
    
    coefficients = cbind(point.est, se.est, conf.lower, conf.upper, p.value)
    
    linear.predictors = va %*% point.est[1:ncol(va)]
    if(param=="RR") param.est = exp(linear.predictors)
    if(param=="RD") param.est = linear.predictors
    if(param=="OR") param.est = expit(linear.predictors)
    
    sol = list(param = param, point.est = point.est, se.est = se.est, cov = cov, 
        conf.lower = conf.lower, conf.upper = conf.upper, p.value = p.value,
        coefficients = coefficients, param.est = param.est, va = va, vb = vb, 
        converged = converged)
    class(sol) = c("brm", "list")
    attr(sol, "hidden") = c("param", "se.est", "cov", "conf.lower", "conf.upper", 
        "p.value","coefficients", "param.est", "va", "vb", "converged")
    
    return(sol)
    
}


## This function is useful for finding limits on the boundary 
## It gives 0.5*exp(x)*(-1+sqrt(1+4exp(-x))) 
## This is bounded between 0 and 1, and takes value (-1+sqrt(5))/2 at x=0 (Some relation to golden ratio) 
## Limits are 0 and 1 as x goes to -infty and +infty respectively 
## The function will never return NaN given a numerical input

getPrbAux = function(x) {
    if ((x < 17) && (x > (-500))) {
        return(0.5 * exp(x) * (-1 + (1 + 4 * exp(-x))^0.5))
    } else {
        if (x < 0) {
            return(0)
        } else {
            return(1)
        }
    }
} 
